﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace WpfApp1
{
     public class CalHeartRate
    {

		double[] m_dSerial=new double[1250];
		double[] m_dClean = new double[1250];
		 int m_nSerialLen;
		//数组求最大值
		public	double Max(double[] pDataIn, int nDataLen)
		{
			double dMax = pDataIn[0];
			for (int i = 0; i < nDataLen; i++)
			{

				dMax = CalMax(dMax, pDataIn[i]);
			}
			return dMax;
		}


		  double CalMax(double x,double y)
		{
			if (x <= y)
			{
				return y;
			}
			else
				return x;
		}
			 
		unsafe public void MovingWindow( double[] pDataIn, double[] pDataOut, int nDataLen, int nWidth)
		{
			for (int i = 0; i<nDataLen; i++)
			{
				pDataOut[i] = 0;
			}

			for (int i = 0; i<nDataLen - nWidth; i++)
			{
				double fSum = 0.0;
				for (int j = i; j<i + nWidth; j++)
				{
					fSum += pDataIn[j];
				}
				pDataOut[i + nWidth / 2] = fSum / nWidth;
			}
		}

		unsafe void Square( double[] pDataIn, double[] pDataOut, int nDataLen)
		{
			for (int i = 0; i<nDataLen; i++)
			{
				pDataOut[i] = pDataIn[i] * pDataIn[i];
			}
		}


		unsafe void Derivative(  double[] pDataIn,  double[] pDataOut, int nDataLen)
		{
			for (int i = 0; i<nDataLen; i++)
			{
				pDataOut[i] = 0;
			}

			for (int i = 4; i<nDataLen - 4; i++)
			{
				pDataOut[i] = (2 * pDataIn[i + 4] + pDataIn[i + 2] - pDataIn[i - 2] - 2 * pDataIn[i - 4]) / 8;
			}
		}

		unsafe void OriginalDerivative(  double* pDataIn, double* pDataOut, int nDataLen)
		{
			for (int i = 0; i < nDataLen; i++)
			{
				pDataOut[i] = pDataIn[i];
			}

			for (int i = 0; i < nDataLen - 1; i++)
			{
				pDataOut[i] = pDataIn[i + 1] - pDataIn[i];
			}
		}


	unsafe double  BandPass10_17(double x)
		{
			double[] b = new double[]{ 0.047991, -0.17255, 0.20164, 0, -0.20164, 0.17255, -0.047991 };
			double[] a = new double[] { 1, -5.15633, 11.3461, -13.6246, 9.41314, -3.54779, 0.570356 };

			double[] in_n = new double[7];
			double[]  out_n=new double[6];

			for (int i = 0; i < 6; i++)
			{
				in_n[i] = in_n[i + 1];
			}
			in_n[6] = x;
			double sum1 = 0;
			for (int i = 0; i < 7; i++)
			{
				sum1 += b[i] * in_n[6 - i];
			}
			double sum2 = 0;
			for (int i = 0; i < 6; i++)
			{
				sum2 += a[i + 1] * out_n[5 - i];
			}
			double sum = (sum1 - sum2) / a[0];
			for (int i = 0; i < 5; i++)
			{
				out_n[i] = out_n[i + 1];
			}
			out_n[5] = sum;

			return sum;
		}

		unsafe public int  GetVCGParameter(double[] dInPut, int nLen, out int  nLastQrsPos)
		{
			for (int i = 0; i < 1250 - nLen; i++)
			{
				m_dSerial[i] = m_dSerial[i + nLen];
			}

			for (int i = 0; i < nLen; i++)
			{
				m_dSerial[1250 - nLen + i] = dInPut[i];
			}
			// 数据长度
			m_nSerialLen += nLen;
			if (m_nSerialLen > 1250)
			{
				m_nSerialLen = 1250;
			}


			// 带通滤波
			double[] dClean;
			dClean = new double[nLen];
			for (int i = 0; i < nLen; i++)
			{
				dClean[i] = BandPass10_17(dInPut[i]);
			}

			for (int i = 0; i < 1250 - nLen; i++)
			{
				m_dClean[i] = m_dClean[i + nLen];
			}

			for (int i = 0; i < nLen; i++)
			{
				m_dClean[1250 - nLen + i] = dClean[i];
			}

			if (m_nSerialLen < 1250)
			{
				for (int i = 0; i < 1250 - m_nSerialLen; i++)
				{
					m_dSerial[i] = m_dSerial[1250 - m_nSerialLen];
					m_dClean[i] = m_dClean[1250 - m_nSerialLen];
				}
			}

			// 微分
			double[] dDeri = new double[1250];  
			Derivative(m_dClean, dDeri, 1250);

			// 平方
			double[] dSqure = new double[1250];
			for (int i = 0; i < 1250; i++)
			{
				dSqure[i] = dDeri[i] * dDeri[i];
			}

			// 滑动平均
			double[] dMove = new double[1250]; ;
			MovingWindow(dSqure, dMove, 1250, 25);

			// 最值
			double dMaxMove = 0;
			for (int i = 0; i < 1250; i++)
			{
				dMaxMove = CalMax(dMaxMove, dMove[i]);
			 
			}

			double dMaxData = m_dSerial[0];
			double dMinData = m_dSerial[0];
			for (int i = 0; i < 1250; i++)
			{
				dMaxData = CalMax(dMaxData, m_dSerial[i]);
				dMinData = m_dSerial[i] >= dMinData ? dMinData: m_dSerial[i];
			}
			List<int> aHeart = new List<int>();
			for (int i = 20; i < 1230; i++)
			{
				if (dMove[i] >= dMaxMove / 6 && dMove[i] > dMove[i - 10])
				{
					double dMax = 0;
					int nPos = i;
					for (int j = i; j < i + 20; j++)
					{
						if (dMove[j] > dMax)
						{
							dMax = dMove[j];
							nPos = j;
						}
					}
					aHeart.Add(nPos);
					i = nPos + 45;
				}
			}


            if (aHeart.Count() > 1 && dMaxData - dMinData > 0.1 && m_nSerialLen > 500)
            {
                double dMeanRR = 0;
                for (int i = 1; i < aHeart.Count(); i++)
                {
                    dMeanRR += (aHeart[i] - aHeart[i - 1] + 1) * 1.0 / (aHeart.Count() - 1);
                }

                nLastQrsPos = 5000 - 1000 * aHeart[aHeart.Count() - 1] / 250;
                return (int)(60.0 * 250.0 / dMeanRR);
            }

            nLastQrsPos = -1;
			return -1;

		}


	  unsafe	 void burg(double[] x, int n, int p, double[] a, out double  v)
		{
			int i, k;
			double r0, sumd, sumn;
			double[] b,  ef,  eb;
			b = new double[p + 1];
			ef = new double[n];
			eb = new double[n];
			a[0] = 1.0;
			r0 = 0.0;
			for (i = 0; i < n; i++)
			{
				r0 += x[i] * x[i] / n;
			}
			v = r0;
			for (i = 1; i < n; i++)
			{
				ef[i] = x[i];
				eb[i - 1] = x[i - 1];
			}
			for (k = 1; k <= p; k++)
			{
				sumn = 0.0;
				sumd = 0.0;
				for (i = k; i < n; i++)
				{
					sumn += ef[i] * eb[i - 1];
					sumd += ef[i] * ef[i] + eb[i - 1] * eb[i - 1];
				}
				a[k] = -2 * sumn / sumd;
				for (i = 1; i < k; i++)
				{
					b[i] = a[i] + a[k] * a[k - i];
				}
				for (i = 1; i < k; i++)
				{
					a[i] = b[i];
				}
				v = (1.0 - a[k] * a[k]) * v;
				for (i = (n - 1); i >= (k + 1); i--)
				{
					ef[i] = ef[i] + a[k] * eb[i - 1];
					eb[i - 1] = eb[i - 2] + a[k] * ef[i - 1];
				}
			}
		 
		}



		unsafe void psd(double[] b, double[] a, int q, int p, double sigma2, double fs, double[] x, double[] freq, int len, int sign)
		{
			int i, k;
			double ar, ai, br, bi, zr, zi, im, re, xre, xim;
			double ang, den, numr, numi, temp;
			for (k = 0; k < len; k++)
			{
				ang = k * 0.5 / (len - 1);
				freq[k] = ang * fs;
				zr =Math. Cos(-8.0 * Math.Atan(1.0) * ang);
				zi = Math.Sin(-8.0 * Math.Atan(1.0) * ang);
				br = 0.0;
				bi = 0.0;
				for (i = q; i > 0; i--)
				{
					re = br;
					im = bi;
					br = (re + b[i]) * zr - im * zi;
					bi = (re + b[i]) * zi + im * zr;
				}
				ar = 0.0;
				ai = 0.0;
				for (i = p; i > 0; i--)
				{
					re = ar;
					im = ai;
					ar = (re + a[i]) * zr - im * zi;
					ai = (re + a[i]) * zi + im * zr;
				}
				br = br + b[0];
				ar = ar + 1.0;
				numr = ar * br + ai * bi;
				numi = ar * bi - ai * br;
				den = ar * ar + ai * ai;
				xre = numr / den;
				xim = numi / den;
				switch (sign)
				{
					case 0:
						{
							x[k] = xre * xre + xim * xim;
							x[k] = sigma2 * x[k] / fs;
							break;
						}
					case 1:
						{
							temp = xre * xre + xim * xim;
							temp = sigma2 * temp / fs;
							if (temp == 0.0)
							{
								temp = 1.0e-20;
							}
							x[k] = 10.0 * Math.Log10(temp);
						}
						break;
				}
			}
		}

		public	List<double> m_vecPSDRR = new List<double>();
		public List<double> m_vecFreqRR = new List<double>();

		public class HRVFFTPara
		{
			public double VLF { set; get; } = 0;//极低频
			public double ALLPSD { set; get; } = 0;//总功率谱

			public double LF { set; get; } = 0;//低频

			public double HF { set; get; } = 0; //高频

			public double LFnormal { set; get; } = 0; //标准化

			public double HFnormal { set; get; } = 0; //标准化

			public double LFHF  { set; get; } = 0; //LF/HF

		}
		public HRVFFTPara hRVFFTPara;

		unsafe  public void  GetPSDRR(List<int> m_vecRR    )
		{

			// 数据清除
			m_vecPSDRR.Clear();
			m_vecFreqRR.Clear();
			hRVFFTPara = new HRVFFTPara();
			// RR个数
			int nCount = m_vecRR.Count();
			if (nCount < 3) return;

			// RR均值
			double dMeanRR = 0;
			for (int i = 0; i < nCount; i++)
			{
				dMeanRR += 1.0 * m_vecRR[i] / nCount;
			}

			double[] pX = new double[nCount];
			bool bAllZero = true;
			for (int i = 0; i < nCount; i++)
			{
				pX[i] = m_vecRR[i];
				pX[i] = pX[i] - dMeanRR;
				if (pX[i] != 0)
				{
					bAllZero = false;
				}
			}
			if (bAllZero)
			{
				pX[0] = 1;
			}
			int nP = nCount - 1<=18 ? nCount - 1:18;// 最小值
			double[] fA = new double[nP + 1];
			double fV=0;
			burg(pX, nCount, nP, fA,  out fV);

			double[] fB = new double[1];
			fB[0] = 1;

			int nFreqNum = 300;
			double[] fX = new double[nFreqNum];
			double[] fFreq = new double[nFreqNum];
			psd(fB, fA, 0, nP, fV, 1.0, fX, fFreq, nFreqNum, 0);

			for (int i = 0; i < (int)nFreqNum; i++)
			{
				m_vecFreqRR.Add(fFreq[i]);
				m_vecPSDRR.Add(fX[i]);
			}

			//计算 功率 0-0.4Hz
			for (int i=0;i< m_vecPSDRR.Count *4/5;i++)
			{
				hRVFFTPara.ALLPSD += m_vecPSDRR[i] * 0.5 / m_vecPSDRR.Count;
			}
			//计算VLF 超低频 0-0.04
			for (int i = 0; i < m_vecPSDRR.Count * 4 / 50; i++)
			{
				hRVFFTPara.VLF += m_vecPSDRR[i] * 0.5 / m_vecPSDRR.Count;
			}
			//计算 LF 超低频 0.04-0.15
			for (int i = m_vecPSDRR.Count * 4 / 50; i < m_vecPSDRR.Count * 3 / 10; i++)
			{
				hRVFFTPara.LF += m_vecPSDRR[i] * 0.5 / m_vecPSDRR.Count;
			}
			//计算 HF 超低频  0.15-0.4
			for (int i = m_vecPSDRR.Count * 3 /10; i < m_vecPSDRR.Count * 4 / 5; i++)
			{
				hRVFFTPara.LF += m_vecPSDRR[i] * 0.5 / m_vecPSDRR.Count;
			}

			hRVFFTPara.LFnormal = hRVFFTPara.LF * 100 / (hRVFFTPara.ALLPSD - hRVFFTPara.VLF);

			hRVFFTPara.HFnormal = hRVFFTPara.HF * 100 / (hRVFFTPara.ALLPSD - hRVFFTPara.VLF);

			hRVFFTPara.LFHF = hRVFFTPara.LF / hRVFFTPara.HF;
		}



		//double[] data = { 25, 8, 15, 5, 6, 10, 10, 3, 1, 20, 7 };
		//int[] index = getPeaksIndex(trendSign(oneDiff(Constant.data)));

	
		 
		
	}
}
